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We evaluate the density of states of the quenched normal modes of ST2 water, 
and their statistical fluctuations, for a range of densities spanning three regimes 
of behavior of a hydrogen bonded liquid: a lower-density regime of random tetra- 
hedral network formation; in the vicinity of a liquid-liquid critical point; and in a 
higher-density regime of fragile glass-forming behavior. For all cases we find that the 
fluctuations around the mean spectral densities obey the predictions of the Gaussian 
orthogonal ensemble of random matrix theory. We also measure the participation 
ratio of the normal modes across the entire frequency range, and find behavior consis- 
tent with the majority of modes being of an extended nature, rather than localized. 

PACS numbers: 63.50.-x, 05.45.Mt, 24.60.Lz, 61.43.Fs 

I. INTRODUCTION 



The nature of the potential energy surface (PES) of liquids and glasses has attracted 
significant attention over the past decade. Numerous computer simulation studies have 
successfully demonstrated that both the thermodynamic and transport behavior of glass- 
forming liquids can be better understood through a consideration of the PES, in particular 
by examining its local minima, or inherent structures (IS), and the nature of the energy 
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A fundamental property describing the topology of the PES at a particular point in config- 
uration space is the matrix of second derivatives (or Hessian matrix) of the potential energy 
with respect to the degrees of freedom of the system. The eigenvalues and eigenvectors of 
the Hessian matrix characterize the normal modes of the configuration under consideration. 



If the configuration is a snapshot of an equilibrium hquid, the modes are called instanta- 
neous normal modes (INM), and the spectrum of eigenvalues will contain both positive and 
negative values. If the configuration is an IS, the modes are called quenched normal modes 
(QNM), in which case all the eigenvalues are positive. Given an ensemble of IS's evaluated 
from equilibrium liquid configurations, the equilibrium probability density for QNM's of a 
given frequency can be found, and corresponds to the vibrational density of states of the 
system in the harmonic approximation. 

While the shape of vibrational spectra of disordered systems in general can vary widely 
from one system to another, commonalities have been discovered using concepts from random 
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ISj. For example, random matrix theory has been used extensively 



for the analysis of statistical fluctuations of the spectra of complex quantum systems such 
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as complex nuclei, quantum chaotic systems and disordered mesoscopic systems 
Study of the spectra of these systems has revealed that the fluctuations around mean spectral 
densities can always be assigned to one of three universality classes, even though mean 
spectral densities themselves may be highly system dependent. 

Notably, in the context of glass-forming systems, recent simulation studies conducted on 
the vibrational spectra of several liquid and amorphous solid systems have shown that the 
fluctuation properties of the vibrational spectra are described by the Gaussian orthogonal 



ensemble (GOE) of random matrices 
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using spectra derived from both INM and QNM properties, although the clearest indications 
seem to have come from using QNM data. All the glass-forming systems tested to date have 
been found to obey GOE statistics 291]. These systems make use of various spherically- 
symmetric atomic interaction potentials and include both bulk liquids and glasses as well 
as clusters, in both two and three dimensions. The results suggest that the statistics of the 
fluctuations of vibrational spectra in all liquids and glasses may obey a common pattern, 
represented by the GOE universality class. 

The identification of universality in the fluctuations of liquid vibrational spectra can assist 
in checking analytical approaches (see e.g. Refs. [30|, 1311) that attempt to evaluate densities 
of state directly from a system Hamiltonian, a fundamental challenge of liquid-state physics. 
The nature of the universality also informs our understanding of the dynamical processes in 
liquids and glasses. If the universality is that of the GOE of random matrix theory, then tjhe 
implication is that the normal modes of the system are extended (rather than localized) 



thus elucidating the nature of the system's motion on the PES. 

In this work, we address the question of whether this proposed universahty extends 
beyond atomic systems to molecular, network-forming liquids. This is an important class 
of glass-forming systems, as it includes water, silica and many other materials that form 
the basis of amorphous solids in a range of contexts, from biology to earth sciences. Here, 
we will employ computer simulations of a model water potential to represent the class of 
network-forming liquids in which tetrahedral coordination dominates at the local level. As 
described below, this will allow us to extend the search for GOE behavior to a range of 
conditions not accessible using simpler atomic potentials. 



II. METHODS 



To model water, we use the ST2 pair potential formulated by Stillinger and Rahman 32|. 
Over the years, there have been numerous studies that have used this model to elucidate the 
complex phenomena found in water, and other liquids having a local tetrahedral molecular 
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38( 1 . and references therein. While other simula- 



tion models exist that more accurately reproduce water properties, ST2 has the advantage 
of displaying a diversity of extreme liquid behavior in one system. With regard to the PES 
of water, the initial study on the IS in water was conducted by Stillinger and Weber jj]. 
Thereafter there have been numerous studies of liquid water dynamics and its relationship 



to the underlying PES using various model potentials 
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We study a periodic ST2 system with a constant number of molecules = 1728 (unless 
otherwise noted), constant volume V, and with temperature T controlled by a Berendsen 

n □ 

thermostat [39]. We conduct conventional molecular dynamics simulations [40| using the 
"leap-frog" algorithm with a 1 fs time step, and where the long range nature of the ST2 
interaction potential is approximated using the reaction field method, with direct inter- 
molecular interations cut off at 0.78 nm. Our molecular dynamics code uses the SHAKE 
algorithm to constrain the motion of the ST2 molecule to that of a rigid body 41|. 

We study three densities, p = 0.83, 0.93 and 1.09 g/cm^ at T = 260 K. These three 
densities are chosen (following the reasoning of Becker, et al. 38|]) to allow us to explore a 
wide range of liquid behavior in ST2 water. A tetrahedral network of hydrogen bonds is 
predominant throughout this density range. At the same time, changes in the hydrogen bond 



network with density are correlated to qualitative changes in thermodynamic and transport 
properties. At the lowest density we study, p = 0.83 g/cm^, the behavior of ST2 water is 
characterized by the emergence of an increasingly defect-free random tetrahedral network 
as T decreases. This is accompanied by the detection of both a density maximum and a 
density minimum 36|], and by a fragile-to-strong dynamical crossover 37|], by following the 
behavior along this isochore. At the intermediate density of p = 0.93 g/cm^, we enter the 
regime at which a liquid-liquid phase separation (between a low density liquid and a high 
density liquid phase) occurs below approximately T = 245 K [33]. At T = 260 K, the 
system is under the influence of critical fluctuations associated with the critical point of this 
phase transition. Finally, at the highest density we study, p = 1.09 g/cm^, the tetrahedral 
structure of the first coordination shell is increasingly disrupted by the penetration of a 
non-hydrogen-bonded fifth neighbor molecule, inducing dynamical behavior more like that 



of a fragile glass-former 



37|. 



In the following, we analyze the vibrational density of states of the QNM's of the system, 
and the statistics of their fluctuations, at each of these densities. The three distinct regimes 
of liquid-state behavior represented by these densities (respectively: random tetrahedral net- 
work, critical fluctuations, and fragile dynamics) provide a range of extreme thermodynamic 
and dynamical conditions in a network-forming molecular liquid for testing the robustness 
of the universal behavior of vibrational spectra that have been reported for simpler atomic 
systems. Our aim is to test if any of these complexities of a tetrahedral network is sufficient 
to induce an exception to this universality. 

To simplify the comparison of our results at different densities, our data is generated 
at T = 260 K. This is low enough to reveal the characteristic behavior represented by the 
three densities chosen. We have also checked that our main conclusions do not change if 
the analysis at each density is carried out at 300 K. It is important to note that we do not 
expect our results to depend on whether or not the liquid is supercooled at these T. We 
choose to study this low T regime simply to ensure that the network-forming character of 
the liquid is prominent. 

We also analyze the normal modes to determine if they are extended (i.e. involve an 
extensive number of particles of the system), or localized to a finite cluster of molecules. 
Previous work has connected extended modes with those obeying GOE statistics, while 
localized ones are thought to be consistent with Poissonian statistics Previous studies of 
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FIG. 1: (a) DOS versus u) of quenched normal modes plotted for various densities at T = 260 K and 
for = 1728. The plot for each density is an average over 100 configurations, and is normalized 
to unit area, (b) Average participation ratios versus u) of the quenched normal modes for various 
densities at T = 260 K and N = 1728. The plot for each density is averaged over 100 configurations, 
(c) Average participation ratios versus w, plotted on a semi-log scale, of the quenched normal modes 
for p = 0.83 g/cm^ and T = 260 K obtained from simulations of two different system sizes, N = 216 
and 1728 molecules. 



network-forming liquids, inc 



are extendec 
ratios [id, 11 



in nature 



12 



13 



14 



18 
15 



uding silica and water, have reported that a majority of modes 
We study this issue here by evaluating the participation 
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regard is consistent with the behavior found for the statistics of the spectral fluctuations. 



III. DENSITY OF STATES 



We evaluate the density of states (DOS) from liquid configurations harvested from our 
simulations. For each density, the system is allowed to equilibrate with the temperature 



controlled by a Berendsen thermostat. The equilibration runs at each T are carried out for 
the time required for the mean squared displacement (MSD) to reach 1 nm^, or for 100 ps, 
whichever is larger. Note that a MSD of 1 nm^ corresponds to an average displacement of 
about 3 times the distance between nearest neighbour molecules. This equilibration criterion 
is the same as that used in Ref. 36| and produces stable and reproducible time series of 
thermodynamic quantities. After equilibrium is attained, the runs are continued, again with 
T controlled using a Berendsen thermostat. From these production runs we generate 100 
distinct configurations, separated in time by 10 ps. We then quench each configuration to 
its corresponding IS using the conjugate gradient method. 

For each IS, we evaluate the Hessian matrix with Bloch wave vector k = 0. The definition 



of the Hessian matrix used here is exactly the same as that given in Appendix B of Ref. [23 1. 
We then diagonalize the Hessian matrix to obtain both the eigenvalues (A) as well as the 
eigenvectors. The eigenvectors correspond to the normal modes of the IS configuration and 
the eigenvalues are related to the frequencies (u) of the obtained normal modes according 
to = V^. The normal modes obtained in this manner are thus QNM's of the liquid. The 
three smallest eigenvalues are zero as a result of translational invariance of the system, and 
the remaining (6A^ — 3) eigenvalues have positive definite values. 

In Fig. 1(a), we plot the DOS function G{uj) for the QNM's, by constructing the his- 
togram of the frequencies uj. To improve the statistics we combine the results from all 100 IS 
configurations available for a particular p and T. To facilitate comparison, the G{u) curves 
for the three densities are normalized so that the area under each is unity. 
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2J]. Our DOS results are consistent with the qualitative 



pattern found in these earlier studies, both for quenched and instantaneous normal modes. 
As in these studies, we find a distinct band of modes below about 400 cm~^ that has been 
shown to be associated with translational degrees of freedom. The modes above 400 cm~^ 
have been assigned to rotational degrees of freedom. We also note that the qualitative form 
of G{uj) varies little across the three densities studied here; this reflects the predominance of 
the tetrahedral network at all three densities. At the same time, given the changes in ther- 
modynamic and dynamical properties with density, the nature of the fluctuations around the 
mean DOS must be separately examined (see next Section) to check for influences connected 
to density change. 



A distinguishing feature of our data for G{uj) is the gap in the DOS at approximately 
850 cm^^. This feature is absent, or only present as a weak shoulder in the data for other po- 
tentials. In a recent study of SPC/E water, this feature was identified as due to a separation 
of the rotational modes associated with each of the principal axes of the water molecule, with 
two axes contributing to the band between 400 and 850 cm~^, and the third contributing to 



that above 850 cm 



25| . Our results suggest that this separation is even more prominent 



in ST2 water, but further study is required to clarify this point. 

Fig. 1(b) shows p{uj), the average participation ratio as a function of u, calculated using 
the eigenvectors for each of the densities. For an individual normal mode a, the definition 
of the participation ratio used here is. 
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where is the vector contribution of the degrees of freedom of particle i to the normalized 
eigenvector e^, that is, = Y^iLi ■ The average participation ratio p{uj) is the average of 
Pa for normal modes with eigenfrequencies within 6uj of a given u, where 600 is the bin width 
used to numerically evaluate p{uj). Defined in this way, p{uj) measures the average fraction 
of molecules in the system involved in modes having a given frequency. Extended modes 
(i.e. modes involving an extensive number of molecules) will have values of the participation 
ratio that are independent of system size, while localized modes will have values that scale 
with the system size as 17l |. 

To test which modes are extended or localized, we present in Fig. 1(c) the participation 
ratios at p = 0.83 g/cm^ and T = 260 K obtained from simulations of two different system 
sizes, N = 216 and 1728 molecules. This data is presented on a semi- log scale in order 
to reveal the behavior of the data at small values of the participation ratio. Since these 
two systems differ in size by a factor of 8, localized modes in the larger system should 
have participation ratios approximately 1/8 of those in the smaller system. We find that 
differences on this scale only occur in the gaps between the bands in the DOS occurring at 
400 and 850 cm~^. Hence, even though the participation ratios near the maximum of the 
highest frequency band are as low as 0.3, these still have behavior consistent with extended 
modes of the system. This analysis illustrates that the magnitude of the participation ratio 
alone can be a misleading indicator of whether or not modes at a particular frequency are 
extended or localized. 



10000 
8000 
6000 
4000 
2000 
0, 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 






^^^Region IV 


Region III/ 
/ 


- 


/Region II 


- 


- /Region 1 

77 






200 400 600 800 1000 1200 
CO (cm'^) 



FIG. 2: The integral H{lo) of the density of states for a single spectrum for p = 0.83 g/cm^ at 
T = 260 K. Also indicated are the four regions within which a quadratic polynomial is fitted to 
this function. 

IV. FLUCTUATIONS 



Our next aim is to examine the statistical fluctuations around the average DOS. Unless 
stated otherwise, the fluctuation properties presented below are computed for the eigen- 
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frequencies, in contrast to earlier studies pj Il5l . 1161 ] that used the eigenvalues. For 
a particular IS, we denote the elements of the obtained spectra of eigenfrequencies by Ua 
with a = 1,2,..., 6A^. Since the present system is periodic, the flrst three elements in the 
spectrum are zero as they are associated with the translations of the entire system along the 
three Cartesian directions, and the remaining (6A^ — 3) positive frequencies are characterized 
by deflning a mean local density as well as fluctuations around it. 

In order to study the fluctuation properties, we flrst must transform the frequencies in 
such a way that the average spacing between two successive frequencies in a given spectrum 



is unity. This process is known as "unfolding" the data 
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unfold the data, we deflne H{u) to be the number of frequencies equal to or less than uj 
(shown in Fig. 2), and S(uj) as a smooth function that passes through if in a best-flt sense. 

nnn 

Note that, as also found in Refs. [13|, IM, 115| in the present case there is no simple function 
that passes smoothly through the whole of H{uj). In order to overcome this complication, 
we divide the complete spectrum into four regions as shown in Fig. 2. Within each region, 
we use a quadratic polynomial D{uj) = a + b u + c uj"^ as an approximation for S{uj). The 
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FIG. 3: (a). Probability density P{s) of normalized nearest neighbor level spacings s for several 
densities at T = 260 K. (b). Variance of the number of levels in intervals of length r for several 
densities at T = 260 K. Symbols are the same as shown in (a). In both (a) and (b), the data are 
compared to the prediction of the GOE (solid line). 



values of a, b, and c are obtained by a standard least-squares fitting procedure. We also 
calculate the misfit (or residual) function 13,[l4, 15] corresponding to the fits in each region 
to check how well D{uj) approximates H{uj). The plots of the misfit functions in our case 
are qualitatively similar to the plot in Fig. 1(b) of [14.]. In order to further improve the 
fit, in each of the four regions we eliminate subregions where the misfit function has a very 
irregular behavior. In the remaining regular subregions we fit a new quadratic function to 
the misfit function and correct D{uj) by these quadratic functions, which yields the desired 
unfolding functions. For each of the fluctuation properties reported here, we combine the 
data from all the fitted subregions of all the spectra for a given density and temperature. 

The first fluctuation property that we report is the distribution P{s) of the normalized 
nearest-neighbor spacings s of the frequencies of the unfolded spectra [Fig. 3(a)]. To improve 
the statistics for each of the three densities studied at T = 260 K, we combine the spacing 
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FIG. 4: (a) Skewness parameter 71 (r), and (b) excess parameter 72(r) of the distribution n(r), 
the number of levels in intervals of length r, shown for several densities at T = 260 K. Symbols in 
(b) are the same as shown in (a). In both cases, the data are compared to the prediction of the 
GOE (solid line). 

data from all four subregions, and from all of the unfolded individual spectra (one for each 
of the 100 IS configurations}. Also shown in Fig. 3(a) is the prediction based on the GOE of 
random matrix theory 26|, 1271]. The agreement of the theory and the data is extremely close. 
In addition, we also check that selection of a random individual spectrum and analysis of 
each of the four regions separately still shows that the fluctuation properties correspond to 
those of the GOE. 

The second fluctuation property we report is S^(r), the variance of the number of levels 
n{r) within an interval of length r located randomly in a given unfolded spectrum. This 
is plotted for r as large as 20 in Fig. 3(b). It must be emphasized that the quadratic 
correction applied to the fitting function D{uj) is very important in calculating S^(r). This 
calculation is extremely sensitive even to very small errors in the aproximation to S. The 
contribution of any such error to S^(r) grows as r^, whereas the GOE prediction for S^(r) 



grows only as ln(r). As described in Ref. 

MM, 

even though the broad contour of the 
misfit fuction is the same for all spectra, the exact locations of the irregular domains vary 
from one spectrum to another and this might be the cause of the observed deviation observed 
towards r = 20. A detailed analysis using a choice of subdomains tailored to each individual 
spectrum IJ, ll5|] would help to clarify this issue. At the same time, we emphasize that the 
agreement with theory shown in Fig. 3(b) is of the same quality as observed in the case of 



amorphous clusters [li 
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In Figs. 4(a) and (b), we plot the skewness parameter. 
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and the excess parameter 
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of the n(r) function. Here M is the number of data values nj sampled from the n(r) functions 
for individual unfolded spectra, a and n are respectively the standard deviation and mean 
of the M values of n(r) sampled. Also included on the same plots are the predictions for the 
GOE. These predictions are calculated on the basis of a large ensemble of 500 x 500 matrices 
belonging to the GOE. We find extremely close agreement with the theoretical prediction 
for each of the densities. 

Finally, we recall the indication from Fig. [T](c) that the modes in the vicinity of 400 
and 800 cm~^ do not involve an extensive number of particles, and are therefore localized 
modes. If so, the spectral fluctuations in these subregions should not conform to the GOE 
predictions. To test this we have evaluated both P{s) and S^(r) using the same procedure 
as used to create the data in Fig. [3], except we instead take the spectral data from the 
region between "Region 11" and "Region III" , and between "Region III" and "Region IV" , 
identified in Fig. [21 The results are shown in Fig. While the data are noisier than in Fig. [3] 
because of the reduced statistics, there is a systematic departure from the GOE prediction. 



While the data do not follow a purely Poissonian form, the peak position o 
the direction expected when Poissonian effects contribute to the statistics 



P(s) is shifted in 



421]. The data are 



thus consistent with the trend toward localization of the modes in these spectral subregions. 
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FIG. 5: (a). Probability density P{s) of normalized nearest neighbor level spacings s for several 
densities at T = 260 K, calculated for the subregions in the vicinity of 400 and 800 cm~^ (b). 
Variance of the number of levels in intervals of length r for several densities at T = 260 K, 
calculated for the subregions near 400 and 800 cm~^. Symbols are the same as shown in (a). In 
both (a) and (b), the data are compared to the prediction of the GOE (solid line). 

V. DISCUSSION 

Previous studies on simple liquids, both for clusters and in bulk, and for two and three 
dimensions, suggest that GOE behavior in the statistics of the DOS is ubiquitous. Our re- 
sults show that this commonality also extends to a molecular liquid across several densities, 
each representative of a distinct regime of complex liquid behavior. In all the regimes we 
examine, neither random tetrahedral network formation, nor proximity to a liquid-liquid 
critical point, nor entering a higher density regime of fragile glass-forming behavior, is suf- 
ficient to induce a departure in the statistics of the DOS from the predictions of the GOE 
of random matrix theory. Our work therefore considerably extends the evidence in support 
of the robustness and universality of this behavior in liquid systems. 



We note that the Hessian matrix we obtain is subject to constraints among its elements 
that have their origin in the conservation of momentum in the associated molecular system. 
Such a Hessian matrix is distinct in character from one not subject to a momentum con- 
servation constraint, such as the Hessian associated with an electronic DOS in a disordered 
system. The differences in the eigenvalues and eigenvectors in these two cases has been 
discussed in Refs. 
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4J| . In particular, a recent analysis indicates that real-symmetric 



random matrices subject to momentum conservation may be fundamentally associated with 
the GOE ensemble [45|. Our results are largely consistent with this analysis. 

Our results also shed light on the nature of localization of the normal modes of disordered 
systems. While the participation ratios found across a significant fraction of the frequency 
domain fall below 0.6, analysis of the scaling with system size suggests that the vast majority 
of the normal modes (excepts those in the "gaps") are extended in nature, rather than 
localized. This is consistent with the fact that these modes obey the statistics of the GOE. 
Hence we have evidence that modes with values of p{uj) as low as 0.3 are still extended 
in nature. This result is consistent with an earlier report on the minimum value of the 
participation ratio for extended modes in a model soft sphere system It would be 

interesting to examine the nature of extended modes that involve such a low fraction of 
particles in the system. We leave this question for future work. 
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